Computation of Eigenvalues and Eigenvectors for the Discrete Ordinate and Matrix Operator Methods in Radiative Transfer

Nakajima and Tanaka showed that the algebraic eigenvalue problem occurring in the discrete ordinate and matrix operator methods can be reduced to finding eigenvalues and eigenvectors of the product of two symmetric matrices, one of which is positive definite. Here, we show that the Cholesky decomposition of this positive definite matrix can be used to convert the eigenvalue problem into one involving a symmetric matrix. The Cholesky decomposition is extremely stable and is expected to improve the speed of the eigenvalue/eigenvector computation. After a careful comparison of the Nakajima and Tanaka procedure, our new Cholesky decomposition method and the original procedure suggested by Stamnes and Swanson, we find (contrary to our expectations) that the Stamnes and Swanson prescription is still the most accurate because it avoids round-off errors due to matrix multiplications needed to symmetrize the matrix in the two other procedures. We also find that, when the QR algorithm (used to solve the asymmetric eigenvalue problem in the Stamnes and Swanson procedure) is changed to avoid complex arithmetic, the speed becomes comparable to that of the two other procedures based on reduction to symmetric matrices.